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Numerical results from a study of boson stars under nonspherical perturbations 



£h ■ using a fully general relativistic 3D code are presented together with the analysis of 

emitted gravitational radiation. We have constructed a simulation code suitable for 
the study of scalar fields in space-times of general symmetry by bringing together 
components for addressing the initial value problem, the full evolution system and the 



detection and analysis of gravitational waves. Within a series of numerical simulations, 
we explicitly extract the Zerilli and Newman-Penrose scalar ^4 gravitational waveforms 
when the stars are subjected to different types of perturbations. Boson star 
systems have rapidly decaying nonradial quasinormal modes and thus the complete 
gravitational waveform could be extracted for all configurations studied. The 
gravitational waves emitted from stable, critical, and unstable boson star configurations 
are analyzed and the numerically observed quasinormal mode frequencies are compared 
with known linear perturbation results. The superposition of the high frequency 
nonspherical modes on the lower frequency spherical modes was observed in the 
metric oscillations when perturbations with radial and nonradial components were 
applied. The collapse of unstable boson stars to black holes was simulated. The 
apparent horizons were observed to be slightly nonspherical when initially detected and 
became spherical as the system evolved. The application of nonradial perturbations 
proportional to spherical harmonics is observed not to affect the collapse time. An 
unstable star subjected to a large perturbation was observed to migrate to a stable 
configuration. 
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1. Introduction 

Spin zero particles play an important role in particle physics and cosmological models. 
Studies of the early universe have suggested that these types of bosons may have played 
significant roles in determining its history. These roles range from the inflaton, a scalar 
field that could be responsible for the inflation of the universe at very early stages pQ, 
up to the quintessence, which is a scalar field used to model the dark energy component 
of the current universe Bosonic particles could come together through some kind of 
a Jeans instability mechanism to form gravitationally bounded objects such as boson 
stars (complex scalar field) or oscillatons (real scalar field) |3j. These objects would be 
held together by the balance between the attractive force of gravity and the dispersive 
effects of the uncertainty principle (the wave character of the scalar field). Boson stars 
in the Newtonian regime have been considered as candidates for dark matter in the 
context of galactic halos jl]. In the strong field regime these stars have been shown to 
be reliable models for the supermassive objects in the center of some galaxies and it 
has been shown that they may play the role of black hole candidates 6j . 

If these objects exist, the process of formation could include multipolar components 
of the gravitational field that should transform into gravitational waves. Due to 
post-formation perturbations, boson stars could generate gravitational radiation. The 
head-on collision of two boson stars has been investigated previously as a source of 
gravitational waves [7j. When two boson stars collide to form a black hole very little 
scalar radiation is emitted and most of the energy is lost in gravitational radiation [Jj. 
Ryan [H] has studied the scenario of a boson star of mass greater than several M & 
coalescing with a particle such as a small black hole. He finds that this system can 
generate a signal detectable by a gravitational wave interferometer such as LIGO and 
future generation gravitational wave detectors. Recently, Kesden et al. jU] have studied 
the infall of a stellar mass compact object onto a supermassive boson star and studied 
an approximate gravitational wave signal, which could point to differences between 
supermassive horizonless objects (like those in Ref. [H]) and supermassive black holes in 
the center of galaxies. 

The detection of boson stars, either as compact objects or as bricks of the dark 
matter component, necessitates studying their behavior under general nonspherical 
perturbations. Extracting their particular gravitational wave signatures, under general 
perturbations, requires a fully general relativistic 3D code. Boson stars also have 
an important role in the field of numerical relativity aside from their relevance to 
astrophysics. The performance of long-term stable 3D simulations of systems with 
general symmetry and construction of a full gravitational wave signal is a timely and 
challenging problem in general relativity. Boson stars are amenable to stable 3D 
simulation due to their smooth surface and the absence of singularities, offering a suitable 
testbed for the field. 

The properties of spherically symmetric boson stars are well studied in the literature 
[TU1 [TTJ [TJJ [T31 EH EE El Ej- Critical solutions and spherically symmetric boson 
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stars brought to the threshold of black hole formation were studied by Hawley and 
Choptuik in Ref. [18J. They performed a comparison between radial quasinormal 
mode frequencies of boson stars calculated via perturbation theory and those of the 
numerically constructed critical solutions. An investigation of the properties of boson 
stars beyond the spherically symmetric case was performed by Yoshida et al. [19J. 
The authors calculate the quasinormal mode frequencies corresponding to nonradial 
pulsations of a boson star. They perform a linear expansion about the spacetime of a 
spherically symmetric static boson star configuration according to the decomposition 
into tensor spherical harmonics developed by Regge and Wheeler. Their results consist 
of a series of complex quasinormal mode frequencies for even parity perturbations for 
three models of stable, critical, and unstable boson star configurations. These modes 
are damped in a short time because of large imaginary parts of the frequencies. 

The focus of this paper is the nonspherical perturbation problem of spherical boson 
stars and the gravitational wave signals generated by such systems. The simulations 
shown here are performed using the fully general 3D simulation code presented in [20J, 
which is based on the Cactus Computational Toolkit [21]. The existence of the rapidly 
damping quasinormal modes predicted by Yoshida et al. is confirmed numerically. Both 
the Zerilli and the Newman- Penrose \l/4 gravitational waveforms are fully extracted for 
the first time for boson stars. The generated waveforms are compared to modefits 
constructed with the £ = 2 quasinormal mode frequencies calculated by Ref. [TH] . After 
the waveforms are extracted numerically, the energy loss to gravitational radiation is 
estimated. 

In addition to subjecting stars to perturbations proportional to £ = 2 spherical 
harmonics that allowed the comparison between our waveforms and the frequency 
spectrum predicted by Ref. P3J, we subjected stars to more physical perturbations 
and studied (1) the collapse of an unstable boson star to a black hole and (2) the 
migration of a star from the unstable to the stable branch. The collapse of an unstable 
branch boson star was accelerated by applying a mixture of nonradial and radial 
perturbations. The apparent horizon was observed to be measurably nonspherical 
when the horizon formed before the complete relaxation of the nonradial modes. The 
migration is initiated by applying a radial perturbation that significantly reduces the 
mass of the star. A nonspherical perturbation is superimposed on the spherical one and 
the system is evolved. At early times (before a single radial oscillation is complete) the 
full gravitational wave signal is extracted. The migration process was further followed to 
track the continuing low frequency radial oscillations of the metric function g rr (lasting 
hundreds of oscillations of the underlying complex scalar field - a very long time scale 
simulation) . 

The paper proceeds by providing a brief overview of the mathematical background 
that describes the evolution system and the formulation of the initial value problem. 
Sec. I3.1l discusses the different types of perturbation methods and the different boson star 
models considered. The results of simulations are then presented in Sec. 13.21 for stable, 
critical and unstable boson star configurations under small nonradial perturbations. The 
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collapse of an unstable star to a black hole and the formation of the apparent horizon are 
subsequently presented in Sec. 13.31 Finally, the migration of an unstable branch boson 
star to the stable branch under explicit radial and nonradial perturbations is detailed 
in Sec. 13.41 The results are summarized in the conclusion. 

2. Mathematical Background 

2.1. General Relativistic Evolution System 

The action describing a self-gravitating complex scalar field in a curved spacetime is 
given by 



vm 2 )]) (i) 



where R is the Ricci scalar, g^ u is the metric of the spacetime, g is the determinant of the 
metric, $ is the scalar field, V its potential of self-interaction, and the geometric units 
G = c = 1 have been used. The variation of this action with respect to the scalar field 
leads to the Klein-Gordon equation for the complex scalar field, which can be written 

as 

dV 

*** ~ = °- (2) 

When the variation of Eq. (1) is made with respect to the metric g^, the Einstein's 
equations G^ u = SnT^ arise, and the resulting stress energy tensor reads 

= \[d,&dy* + d^d^*] - \gA^^« + m$| 2 ))]- (3) 

In the present manuscript we focus on the free-field case, for which the potential is 
V = m 2 |$| 2 , where m is interpreted as the mass of the field. 

In order to find solutions to the Einstein-Klein-Gordon system of equations we use 
the 3+1 decomposition of Einstein's equations, for which the line element can be written 

as 

ds 2 = -a 2 dt 2 + (dx i + Fdt) (dx j + (3 j dt) (4) 

where 7^ is the 3-dimensional metric; from now on latin indices label the three spatial 
coordinates. The functions a and (3 1 in Eq. (J3J are freely specifiable gauge parameters, 
known as the lapse function and the shift vector respectively. The determinant of the 3- 
metric 7 is defined as 7 = 7y7 y . Throughout this paper the standard general relativity 
notation is used. The Greek indices run from to 3 and the Latin indices run from 1 
to 3. 

The Klein-Gordon equation can be written as a first-order evolution system by first 
splitting the scalar field into its real and imaginary parts: $ = <pi + i<p2, and then 
defining eight new variables in terms of combinations of their derivatives: II = tti + ivr 2 
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and if) a = ipi a + iip2a with n i = (^y/a)(d t (f)i — f3 c d c (pi) and if)\ a = d a (pi and similarly 
(1 — > 2). With this notation the evolution equations become 

(X 

dt<f>i = — 7ri + /^lj (5) 

7 2 

a 

^tVla = d a (— 7T X + Wlj) 
72 

1 dV 
dt-Kx = dj(ayfi$) - 2 a VT^p0i 

and (1 — > 2). On the other hand, the geometry of the spacetime is evolved using 
the BSSN formulation of the 3+1 decomposition. According to this formulation, 
the variables to be evolved are \& = Iii(7y7 v ')/12, 7^ = e _4 *7y, = -f^K^, 
Aij = e~ m (Kij — jijK/3) and the contracted Christoffel symbols P = ^ k Y l - k , instead 
of the usual ADM variables 7^ and Kij. The evolution equations of these new variables 
are described in Refs. 

dtV = - \uK (6) 


dtjij = -2aA ij (7) 



d t K = -Y j DiD ja + a 



+ ^K 2 + l - {-T\ + T) 



(8) 



dtAj = [-DiD ja + a {R %J - T %J )} TF (9) 
+ a {KA l3 - 2A U A 1 ^ (10) 



V = -2A ij a j + 2a(f i jk A kj 



^ 27^^ + 1 W,). (11) 



<9^' V ' >' ' ^ ' m 3 
where -Dj is the covariant derivative in the spatial hypersurface, T is the trace of the 
stress-energy tensor (JHJ) and the label TF indicates the trace-free part of the quantity in 
brackets. The coupling between the evolution of the BSSN variables and the variables 
describing the evolution of the scalar field is first order. That is, Eqs. (5) are solved using 
the method of lines with a modified version of the second order iterative Crank-Nicholson 
(ICN) integrator (see Ref. [23]). After a full time step the stress-energy tensor in Eq. 
(3) is calculated and used to solve the BSSN evolution equations with an independent 
evolution loop based on the standard second order ICN [23] . 

In addition to the formulation of Einstein's equations and the first order form of 
the Klein-Gordon equation we use certain gauge choices to determine the lapse function 
of Eq. (4) (the shift vector is zero in all of our simulations). The evolutions presented 
in Sec. 13.21 were carried out using the 1 + log slicing condition with the lapse given by 
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d t a = —2aK, where K is the trace of the extrinsic curvature. The more dynamic 
evolutions in Sees. 13.31 and 13.41 black hole formation from boson star collapse and 
migration of an unstable boson star to the stable branch, used a combination of both 
maximal and 1 + log slicing. Maximal slicing was enforced using a K-driver j2E] • 

In order to set up the correct scaled quantities to be evolved we use dimensionless 
variables. For an equilibrium configuration the complex scalar field has a time 
dependence of the form $(r, t) = 0(r)e^°*. This implies that the stress energy 
components given in Eq. (3) and therefore the geometry of the spacetime have to be 
time-independent, whereas the field is oscillating with the characteristic frequency u 
that depends on the central value of the field 0(0). This characteristic frequency ujq 
together with the mass of the boson m provide the natural dimensionless units, 



^code = mr/M^, t code = u t, (12) 



/47T m 

'code t, r T*> C^code 



M P1 r ' " coae M> ' 

which are the ones used within the present analysis. For further code details refer to 
Ref. [20]. 

2.2. Formulation of the Initial Value Problem (IVP) 

Performing numerical simulations in the 3+1 ADM formulation requires the construction 
of valid initial data on some constant time hypersurface. This task entails finding 
a numerical solution to the four nonlinear coupled elliptic equations given by the 
Hamiltonian and momentum constraint equations of the ADM formalism. Finding 
solutions for physically interesting systems that are not highly symmetric can be very 
difficult and computationally expensive. The presence of matter terms as a source 
increases the complexity even further. 

We consider the case of time symmetric initial data where the extrinsic curvature 
Kij = on the initial hypersurface (t = 0) and a background scalar field configuration 
<f>(r) is assumed. Under the condition of time symmetry the equations decouple and 
it becomes sufficient to solve only one nonlinear elliptic equation, i.e., the Hamiltonian 
constraint. For a massive complex scalar field the constraint equation takes the form 

fl=167rp(7 O -,0,n) (13) 

where the scalar curvature is R = j^Rij and p is the energy density. Using the York 
formalism, we introduce a conformal factor \& that relates the actual 3-metric 7^ to an 
initial guess 7^- through the relation 7^ = ^7^. Under this transformation the scalar 
curvature is given by 

R = l-R-1-A^ (14) 

In the standard York procedure, all the terms within the energy density are transformed 
uniformly using the relation p = ^~ 8 p. This transformation would require finding a 
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new bosonic field configuration corresponding to the new energy density p, which is a 
nontrivial problem. We circumvented this difficulty by explicitly writing the dependence 
of the energy density on the metric functions. The conformal transformation is applied 
term by term within p, leaving the scalar field configuration unchanged. On the initial 
hypersurface the Hamiltonian constraint (Eq. (JEJ)) can be rewritten as the nonlinear 
elliptic equation 



\]>4 



\]>5 



16tt 



1 1 



1 1 



2* 12 7 1 2^ 4 1 



--m 2 <t>\ + 



(15) 



We start with an initial guess and then solve the constraint equation for the conformal 
factor This formulation of the IVP problem was first used in Balakrishna's Ph.D. 
thesis for a variety of physical situations related to boson stars |7j. 



3. Evolution of Perturbed Boson Stars 



Table ^ details the four different boson star configurations studied in this paper. The 
mass profile plot in Fig. ^ shows the mass versus central field density (f)(0). All 
configurations listed in Table ^ have negative binding energy (M < N m). Model 1 
is a star on the stable branch, Model 2 is a critical boson star, and Models 3a and 3b 
both lie on the unstable branch. These four models are perturbed in different ways, 
evolved, and analyzed as described below. 





Boson St 


ar Configurations 


Model 


0(0) 


M 


-R95 


1 


0.1414 


0.584 


0.91 9.10 


2 


0.271 


0.633 


0.85 6.03 


3a 


0.3536 


0.622 


0.82 5.05 


3b 


0.300 


0.631 


0.84 5.65 



Table 1. Physical characteristics of the boson star configurations studied in this paper 
are listed. The field 4> is i n units of (l/-\/47r)Mpi- The unperturbed mass M of the 
boson star is in units of Mpj/m. The equilibrium scalar field frequency uj is in units 
of m/Mpj and the radius R95 has units of Mpj/m. By the radius of the boson star we 
are referring to the 95% mass radius of the star. Mpi and m are the Planck mass and 
the mass of the scalar field, respectively. 



All simulations presented in this paper are performed using octant symmetry. More 
general simulations on a full grid are the subject of future work. This will allow the 
study of the effect of using octant symmetry on the accuracy and stability of boson 
star simulations. The waveforms are calculated using the Cactus Extract module for 
the Zerilli function and the Cactus Psikadelia module for ^4 with a radial tetrad 
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choice [21j. This particular wave extraction code has been widely used in the past 

E3 123 123 Em ED 132 1331 ■ 

The Newman-Penrose complex scalars \I/o ^4 represent the ten independent 

components of the Riemann tensor projected onto a basis formed by a null tetrad. These 
scalars have asymptotic behavior \I/ n ~ r ™-5_ Thus, the scalar ^4 dominates in the 
distant wave zone and approximates an outgoing gravitational wave. 

The Zerilli function j3El is a gravitational waveform that is generated in the study 
of linear perturbations of the Einstein field equations. When the full system of the 
Einstein equations is expanded about a background Schwarzschild metric with a small 
perturbation function (proportional to spherical harmonics Yi m ), the resulting linear 
equations contain a set describing even parity perturbations. For the Zerilli function 
ipz, the energy loss to gravitational waves can be calculated using 

Mass Profile and Selected Models 

0.7 




Figure 1. The mass versus central field density cf>(0) for ground state boson star 
configurations is shown. The mass increases to a maximum of 0.633Afpj/rn. This is 
the critical boson star configuration that separates the stable and unstable branches. 
The points shown in the curve correspond to the particular equilibrium configurations 
analyzed below. 



3.1. Methods of Perturbation 

Various initial data sets that represent different classes of nonspherical perturbations of 
boson star configurations are prepared. The equilibrium field configuration, spherically 
symmetric solution $(r), is taken as a base and different types of perturbations are 
added to this configuration. One method involved a scalar field perturbation of the 
form 

<5$ = £t m Y em (6, ¥>)/(r)$(r) (17) 

im 
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x 



Figure 2. Conformal Factor isosurface W = 1.00001 after the initial value problem 
solve. The perturbation of the field is proportional to the Y20 spherical harmonic. The 
dumb-bell shape of the harmonic is apparent in the shape of the conformal factor and 
thus in the initial corrections to the metric. 



where f(r) is a weighting function. 

Various weighting functions were studied. 

• Uniformly weighted: For this perturbation the weighting function was constant, 
/(r) = 1. In this case the scalar field was perturbed across the width of the star, 
and was of significant size near the origin where the field density cf) is maximum. 

• Parabolic weighted : 



(where R p is a radius at which the perturbation is centered). In this case the 
perturbation is set to zero at the origin. This results in a more slight perturbation, 
away from the central peak of the field density of the star. 

More general perturbations involved both radial and nonradial components. One 
such radial-nonradial (R-NR) perturbation type was constructed by giving full weight 
to a spherically perturbed profile along one direction, while keeping an unperturbed 
profile in the perpendicular directions and a smooth linear combination of the two in 
the intermediate region. 

For all the perturbations described in this paper the metric functions are not 
explicitly perturbed. The initial value problem solver starts with the perturbed scalar 
field configuration described above and an initial guess for the metric functions and then 
solves the constraint Eq. (fl3j) for the corrected metric functions. The initial guess for the 
metric functions is taken to be the usual spherically symmetric models of Ref. P ^ ITI H ITT] 
obtained by solving a one-dimensional eigenvalue problem for time independent metric 
functions and for a field configuration oscillating with fixed frequency u and constant 
amplitude at the origin. These perturbations are similar to those performed by Ref. 
[T3j . Alternative perturbations could involve external fields that evolve independently 
of the scalar field of the star (see Ref. [IB])- 




(18) 
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Amplitude of VP 4 vs. Resolution 
Model I Stable Star - Uniformly Weighted Perturbation 



2.7 




1.7 r 1 



/ &x=0.48 
» Ax=037S 



1.6 



0.2 0.4 0.6 0.8 



1.5 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Resolution Ax 



1.1 1.2 1.3 



Figure 3. The amplitude of the Newman-Penrose scalar ^ is displayed as a function 
of resolution for a collection of simulations of a Model I stable star. In each simulation 
the waveform was extracted at a detector located at r = 56.8. The $4 amplitude 
of a given simulation is represented by black dot on the plot. The inset shows the 
first derivative of the amplitude curve. It can be seen that for a resolution better 
than Ax = Ay = Az = 0.48 the curve flattens and its derivative goes sharply 
towards zero. A linear extrapolation towards infinite resolution shows the amplitude 
at Ax = Ay = Az — 0.375 to vary from the continuum limit by only about 1-2%. 



An example of the results obtained from the IVP solver is shown in Fig. An 
isosurface of the conformal factor for a perturbation of the field proportional to the I20 
spherical harmonic is visualized. The prominent dumb-bell shape of the harmonic along 
the z-axis resurfaces in this display of the corrections to the initial metric. 

3.2. Perturbation cases studied 

3.2.1. 'Pure' nonspherical perturbations The evolution of boson stars under 'pure' 
nonspherical perturbations (5M(r) ~ where 5M is the mass change) that are 
proportional to spherical harmonics was investigated. Under such perturbations 
the quasinormal modes corresponding to nonradial oscillations were found to have 
frequencies with both real and imaginary parts [19J (even for configurations that are 
unstable or critical in nature with respect to spherical perturbations). Thus, it was 
anticipated that equilibrium boson star configurations would be stable with respect to 
nonspherical perturbations. 

A series of simulations of perturbed star configurations was performed by taking 
different values for the coefficients e^ m and the weighting function fir) in Eq. (fTTj) . 
These perturbations do not significantly change the mass of the star or the number 
of particles (even for a large nonspherical perturbation the mass change is less than 
0.01%). The gravitational wave content is calculated using the even parity quadrupole 
(£ = 2) Zerilli waveforms and the Newman-Penrose scalar ^4. The waveforms are 
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Figure 4. A stable boson star (Model 1, Table is subjected to a small uniformly 
weighted perturbation proportional to Y20 with £20 — 0.032. (a) The Zcrilli £ = 2, m — 
and (b) the '5 4 waveforms are shown in agreement with modcfits (for time intervals 
[65.0, 79.2] and [71.6, 81.9], respectively) using the first twelve quasinormal modes of 
Ref. |19|. Both waveforms exhibit a low frequency precursor that likely represents the 
initial shaking of the star due to an instantaneous finite perturbation. ^4 appears to 
be a more sensitive indicator of the initial dynamics of the star, (c) Coefficients of the 
modes for the two waveforms are shown to have similar profiles. The contribution of 
the first two modes is negligible, probably due to the failure of the WKB approximation 
in determining these modes accurately |19j . 



computed at several different detector locations in the spacetime. The detectors had 
to be located sufficiently far from the center of the star for the ^4 scalar to dominate 
the other scalars in the Newman- Penrose formalism. Since the boson field tails off 
exponentially, the spacetime in the exterior regions (r > i? 95 ) becomes increasingly close 
to Schwarzschild as the radius increases. The perturbation formalism associated with 
the Zerilli waveforms is valid when the metric functions describing the spacetime are only 
small perturbations away from those of the Schwarzschild spacetime. Thus detectors 
measuring Zerilli waveforms should be located at a distance greater than several -R95 
from the center of the star. Simultaneously, the detectors cannot be too close to the 
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boundary in order to avoid the effects of reflection from the edge of the grid. These two 
requirements coupled with computational limitations on the grid size and the need for 
good resolution are the main challenges of our numerical simulations. 

A uniformly weighted perturbation was first considered. Two different 
perturbations were applied to a Model 1 (Table QJ) stable configuration. One 
perturbation was proportional to the Y20 spherical harmonic with constant of 
proportionality €20 = 0.032. The second perturbation was a linear combination of 
Y20 an d 5^22 with coefficients of €20 = 0.032 and £22 = 0.026. In the beginning several 
simulations were performed to determine the scale of accurate grid resolution. Fig. |3] 
shows the amplitude of the real part of the Newman-Penrose ^4 scalar for the same 
grid size at different resolutions for the stable star configuration under the perturbation 
proportional to Y20 described above. The waveform is extracted at a radial location 
of r = 56.8. From the figure it can be seen that at coarse resolution there is a strong 
dependence on amplitude. However, for resolutions better than Ax = Ay = Az = 0.48 
the curve flattens, with the slope heading rapidly towards zero. The amplitude at the 
best resolution Ax = Ay = Az = 0.375 varies from an extrapolation to the infinite 
resolution limit by only about 1-2%. The final simulation was carried out on a 164 3 
grid with a resolution of Ax = Ay = Az = 0.375. The L2-norm of the Hamiltonian 
constraint for this run remained below 1.2 x 10 -4 for the duration of the evolution. The 
L2-norm of the momentum constraints is about an order of magnitude lower than that 
of the Hamiltonian constraint with a maximum value of 4 x 10 -5 . 

Figs. EJa) and (b) show the £ = 2, m = Zerilli function and the Newman-Penrose 
scalar \l/ 4 from a detector at r = 56.8. The plots also show a linear least squares 
regression using the first twelve I = 2 quasinormal modes calculated in Ref. ^HJ. It 
is observed that after an initial precursor, the star rings into a linear combination of 
its quasinormal modes. The waveforms are observed to be four orders of magnitude 
above the level of the noise, which is taken to be the signal obtained from a simulation 
of an unperturbed spherical star configuration with the same grid size and resolution. 
Although the perturbation of the field contains only i = 2, m = spherical harmonics, 
it is observed that higher order i = 4 Zerilli modes of the gravitational radiation are 
also triggered. However, these Zerilli modes have sufficiently low amplitude (about an 
order of magnitude smaller than I = 2 modes) that they do not contribute significantly 
to the energy loss in gravitational radiation. The waveform ^4 contains all nontrivial 
modes of the star. Nevertheless, the 1 = 1 modes dominate and we are able to fit the 
signal with the I = 2 frequencies determined by Ref. [19] . Using Eq. ()16|) the energy 
loss in gravitational waves was calculated to be about 0.1% of the mass of the star. 

In Fig. |3^c) the relative weights of the twelve modes within the Zerilli and ^4 
waveforms for the stable star from Fig. Efa) and (b) are shown. Both waveforms have 
roughly the same mode content. The slow increase of the imaginary part of the modes 
as the order of the mode increases necessitates the use of numerous modes in our fit. 
Note that in black hole systems, it is usually sufficient to fit waveforms dominated by 
quasinormal modes using only the lowest two modes [37] . In this case, however, it is 
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Figure 5. A Model 1 star is perturbed with a uniformly weighted perturbation 
proportional to Y20 with 620 = 0.032. Although no explicit radial perturbation was 
applied, numerical perturbations due to grid discretization (resolution Ax = Ay = 
Az = 0.375) impose a slight radial perturbation on the system. The amplitude of 
this oscillation is seen to converge toward zero as the resolution is improved, (a) 
The maximum value of g rr is plotted against time. As a result of the nonsphcrical 
perturbation the metric exhibits a high frequency oscillation superimposed on the low 
frequency one. By t = 80 the star has undergone only half of its radial oscillation, while 
the nonradial oscillations have damped off. (b) The difference between the density p in 
the x and z directions is shown as a function of time. This difference is an indicator of 
the degree of asymmetry in the system. It is clearly decreasing as the system evolves 
and the star becomes spherical after the relaxation of the nonradial modes. 



justified to neglect only the modes beyond about twelve, on the grounds that by then 
the imaginary part is sufficiently large and they damp out quickly. As can be seen in 
the figure, even the 11 th and 12 th modes hardly contribute. It is noteworthy that the 
first two quasinormal modes are almost absent from the fit. However, these first modes 
have been described as inaccurate in Ref. JH] by the authors because of the break down 
of their WKB approximation in that domain. It should also be noted that values of the 
eigenfrequencies of Ref. JIH] are dependent on their choice for a definition of the surface 
of the boson star. There is no unique definition for this radius, since the bosonic field 
extends to infinity with exponential damping. Hence, there is a degree of arbitrariness 
in the position of the surface and in the exact numerical values of the eigenfrequencies. 
In their work the surface of the star was artificially placed where the value of the field 
was of order 10~ 5 . We take more realistic field configurations, which are limited by the 
size of the grid and the resolution of our simulations. 

The maximum of the radial metric function g rr is shown in Fig.E^a). The nonradial 
quasinormal mode oscillations are superimposed on the radial oscillation. The latter is 
introduced by grid discretization and its amplitude is resolution dependent. In the 
continuum limit one should observe only the high frequency oscillations. Ref. [20J has 
shown that the amplitude of the radial oscillation converges toward zero as the resolution 
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Figure 6. (a) Figure displaying the t — 2 Zerilli waveforms for a Model 1 stable 
star under a small uniformly weighted perturbation proportional to both the Y20 and 
Y22 spherical harmonics with amplitudes £20 = 0.032 and £22 = 0.026. The frequency 
of the Zerilli I — 2, m = 2 is identical to the frequency of the I — 2, to = Zerilli. 
(b) The ^4 waveform for the same star for the perturbation proportional to a linear 
combination of Y20 and Y22 is compared to the waveform for a Y20 perturbation with 
£20 = 0.032. The very similar frequencies indicate that this configuration has specific 
nonradial quasinormal mode signatures. It can be seen that the "J 4 signal for the mixed 
perturbation is larger because the Newman-Penrose scalar incorporates more modes. 



is improved and that in the spherical case no high frequency oscillations appear. For 
this run at a resolution of Ax = Ay = Az = 0.375, the maximum of the metric has 
risen only from about g rr m&x = 1-166 to g r rmax = 1-170. The frequency of a radial 
oscillation of this star configuration was calculated for this simulation to be about 0.010 
m/Mp[, which agrees with the result ~ 0.0097 m/Mpj obtained from a ID code for the 
same configuration. Since the radial oscillation has a large time period, only half of this 
oscillation has occurred by the end of this run. The damping of the nonradial modes 
can be seen towards the end of the simulation, indicating that the star was settling 
into a spherical configuration. Another simulation that doubled the amplitude of the 
nonspherical perturbation €20 was performed. It was observed that the amplitude of the 
nonradial oscillations in the metric g rr also doubled. This is consistent with the result 
of Ref. ^Hj that the metric perturbation is comparable to the scalar field perturbation. 

Fig. Efb) shows the difference between the density p in the x and z directions as a 
function of time. This difference is seen to be decreasing by an order of magnitude in a 
short timescale. Furthermore, the radial location of the peak of the difference function 
is observed to be moving out as the star becomes more spherical, coinciding with the 
emission of gravitational radiation. 

We now apply a different perturbation to Model 1 and study its evolution. In 
this case a uniformly weighted perturbation proportional to a linear combination of Y20 
and Y22 (amplitudes €20 = 0.032 and €22 = 0.026, respectively) was applied. Fig. EJa) 
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Figure 7. The Zerilli waveform for an unstable boson star (Model 3a, Table under 
a small uniformly weighted perturbation proportional to the Y20 spherical harmonic 
(£20 = 0.032) is displayed. The waveform is shown along with modefit applied over the 
interval [49.2, 60.5]. 



presents the extracted Zerilli i = 2 (both m = and m = 2) waveforms for a stable 
star (Model 1, Table Q). These waveforms are observed to have the same frequency. 
Furthermore, the Zerilli i = 2,m = signal for the star under this mixed Y20 and Y22 
perturbation has the same size and frequency as the t = 2, m = Zerilli function of 
the same star under a perturbation proportional only to Y 2 q (in both cases €20 = 0.032). 
This result shows that the addition of other harmonics to the perturbation does not 
alter the original signal and demonstrates the accuracy of the code. Fig. EJb) shows a 
comparison between the Newman-Penrose scalars ^4 for the pure Y20 perturbation and 
the mixed ^20^^22 perturbation. The signals have the same frequencies, but the signal 
for the mixed perturbation is larger in size because it incorporates more modes. 

Perturbations proportional to the Y20 spherical harmonic with amplitude €20 = 
0.032 were studied in the case of a critical boson star (Model 2, TableHJ) and an unstable 
boson star (Model 3a, Table [TJ. Since these stars are more compact than the stable 
star, they need better resolution. However, they have a smaller radius and do not need 
the same grid size. A 164 3 grid with a resolution of Aa; = Ay = Az = 0.25 was used 
for both these simulations. For the time interval relevant to our study (that prior to 
and encompassing waveform extraction) the simulations proceed accurately with the 
L2-norm of the Hamiltonian constraint remaining below 4.6 x 10~ 4 for the critical case 
and below 6.1 x 10~ 4 for the unstable case, respectively. At late times both stars collapse 
to black holes due to radial perturbations introduced by grid discretization errors. 
More detailed simulations of unstable star collapse with apparent horizon analysis are 
presented in Sec. 13.31 Fig. [7| displays the t = 2, m = Zerilli waveform and modefit for 
this unstable star configuration with a detector located at r = 35. After a fairly large 
precursor, the star settles into a linear combination of its quasinormal modes with late 
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Figure 8. Newman-Penrose scalar ^4 for a Model 1 star under a R-NR perturbation 
is superimposed on the signal due to the Y20 perturbation shown in Fig. E£b) ■ The 
amplitude of the latter has been magnified by a factor of 100 to facilitate the 
comparison. The R-NR perturbation produces a much larger signal. The star did 
not ring in its quasinormal modes initially, but appears to move into these modes after 
a few oscillations. The rapid damping of the gravitational wave signal characteristic 
to boson stars can be clearly seen. 



time agreement between the modefit and numerical waveforms. The extended precursors 
are representative of complex dynamics within unstable systems. 

3.2.2. Perturbation weighted along a chosen axis Perturbations with radial and 
nonradial components as described in Sec. Kill were also studied. In Ref. |Zj it was 
shown that for this type of perturbation the dominant radiation is in the form of scalar 
radiation. The scalar radiation emission process, denoted gravitational cooling [38J, 
relaxes the system to an equilibrium configuration. In this case the Zerilli waveform 
appeared noisy and did not exhibit the quasinormal mode oscillation. The Zerilli 
perturbation expansion about a background Schwarzschild spacetime was probably 
affected when scalar radiation moved through the region of the detector. Even when the 
initial configuration has the detector in an almost Schwarzschild region, the dynamics of 
the star can push scalar radiation into this region during the evolution. The Newman- 
Penrose scalar was observed to be more robust in the presence of scalar radiation. 

Fig. |S1 shows the \l/4 waveform for this perturbation for a Model 1 boson star 
superimposed on the waveform for the same star under the Y20 perturbation of Fig. 
0(b). The R-NR simulation was conducted on a 164 3 grid with resolution Ax = Ay = 
Az = 0.35. Since the Y20 perturbation is significantly smaller, it gives a much smaller 
signal. In order to make a clear comparison the waveform for the Y20 perturbation has 
been magnified by a factor of 100. Although the R-NR perturbation is not huge, it 
is large enough to prevent the star from ringing into its quasinormal modes initially. 
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However, it appears to move into these modes after a few oscillations as can be seen 
in the figure. The oscillation is very clearly damping out on a short timescale as it is 
characteristic for boson stars [TH] . 

3.3. Unstable Star Collapse to a Black Hole 

Horizon Mass 

Unstable Star - Model 3b - Collapse to Black Hole 

0.66 

0.62 
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Time After AH formation 

Figure 9. The mass of the apparent horizon vs. the time elapsed after horizon 
formation is displayed for star collapse simulations Run 1 and Run 2 described in the 
text. Run 1 has an amplification factor of 1.01 and an initial mass of M — 0.642Afp 1 / m. 
Run 2 has radial amplification factor of 1.03 and initial mass M — 0.665Afp[/m. In 
each case the horizon mass is observed to steadily increase towards the total mass of 
the initial configuration as more scalar matter falls through the horizon surface. An 
asymptotic approach to the value of the initial mass is expected and is observed at 
early times, but large errors due to grid stretching cause the mass to rise quickly at 
late time. 

In this section simulations that follow the collapse of unstable boson stars under 
nonspherical perturbations to the formation of black hole horizons are presented. The 
unstable star designated as Model 3b in Table [T] is used as the starting configuration. 
A parabolic weighted nonradial perturbation centered at R p = 2.5 (about half of 
the 95% radius of the unstable star) was applied. This nonradial perturbation was 
superimposed on different radial perturbations that accreted additional matter onto the 
star by means of a uniform amplification of the field of the form /U0 eq (where </> eq (?") is 
the initial equilibrium field configuration). The addition of scalar matter by a spherical 
amplification hastens the collapse process. In contrast, the star is observed to collapse 
at approximately the same rate with or without the nonradial perturbation. 

The time dependence of the mass of the apparent horizon for two collapse 
simulations is displayed in Fig. El The nonspherical perturbations for the systems are 
parametrized by €20 = 0.128 and 622 = 0.104 (using Eq. (fTTj) ). The first simulation, 
denoted Run 1, has a radial amplification of fi = 1.01 while the second, denoted Run 2, 
has a radial amplification of /1 = 1.03. For these amplifications the mass of the star 
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Figure 10. The (£ = 2, m = 0) Zerilli waveform for the case of the unstable star of 
Model 3b with amplification factor of /i = 1.03 and nonradial perturbation specified 
by e 2 o = 0.128 and e 22 = 0.104. 

changes to M = 0.642Mpj/m for the /x = 1.01 case (an increase of about 1.7% from 
the unperturbed mass) and M = 0.665Mp[/m for the \i = 1.03 case (an increase of 
about 5.4% from the unperturbed mass). The horizon is observed to form at a time of 
t = 38.3 within Run 2. This is earlier than the formation time of t = 54.7 in Run 1, 
which is appropriate because of the larger initial mass. It can be seen that the horizon 
mass is 0.463M|[/m when first detected for the case of Run 2. The evolution is carried 
forward for an additional time interval t — 2.1 after horizon formation, corresponding 
to approximately t = 3M of evolution of the black hole. A smooth asymptotic approach 
towards the value of the initial mass is observed at early times after the formation. 
However, at late times the errors in the simulation are becoming large and this is 
manifested in a slight upward turn of the curve along a path that would apparently 
cross the axis corresponding to the initial mass. These features are also echoed in the 
results for Run 1, where the mass of the apparent horizon takes the value 0.427Mp 1 /m 
at formation and increases in an analogous manner. The late time errors are caused by 
the familiar problem of grid stretching associated with singularity avoiding time slicing 
[3^1 HOI HI]. It should be noted that these simulations did not make use of advanced 
shift conditions or excision techniques. Future work will include the employment of such 
methods to extend the duration of the simulations for long times after the black hole 
formation. 

The geometry of the apparent horizons formed in Run 1 and Run 2 is investigated 
by calculating the ratio of the polar to equatorial circumference C r = C p /C e for two 
different great circles ip = and (p = tt/2. The values for C r ((p = 0) and C r (<p = 7r/2) 
are very close to unity, with values in the range C r = 0.9998 ± 0.0002 on all timesteps 
for which horizons are detected in the simulations. The nearly spherical shape of the 
horizons can be understood by considering the Zerilli waveform of Run 2 displayed 
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in Fig. El The waveform starts to ring down near t = 30 and damps to a small 
value by t = 38, which precedes the formation time of the apparent horizon within the 
simulation. The nonradial perturbations have decayed away and the system has become 
nearly spherical by the time of horizon formation. However, it can be demonstrated that 
a slightly nonspherical horizon may be obtained when the formation occurs before the 
full relaxation of the nonradial modes. Another simulation, Run 3, takes an initial 
configuration with a very large amplification factor fi = 1.15 and the same nonradial 
coefficients (€20 = 0.128, €22 = 0.104) as Runs 1 and 2. The initial mass of this 
configuration is 0.801M|j/m, and the subsequent evolution results in a horizon formed by 
t = 20.5. From the values of the circumference ratios displayed in Table El it is seen that 



AH circumference ratios 


^formation 


^cnd 


C r (ip = 0) 


= 1.0017 C r (ip = 0) = 1.0008 


C r {<p = f ) 


= 1.0009 C r (tp = f ) = 1-0004 



Table 2. The ratios of the polar to the equatorial circumference for two different great 
circles (ip — and ip = n/2) are displayed at the horizon formation time and at the 
end of Run 3. The termination time is t cn d — iformation + 3.2 m iformation + 4M, where 
the apparent horizon forms at the time if or mation = 20.5. 

the horizons are measurably nonspherical. The circumference ratios are seen to steadily 
decrease, suggesting the eventual evolution towards spherical symmetry. Furthermore, 
we see that the perturbation proportional to Y22 has resulted in a dependence in the 
ratio C r on the angle ip at which the polar circumference is measured. The formation 
and study of horizons with complex shape from boson star collapse will be investigated 
further by using advanced methods for the evolution of black holes. 

3. 4- Migration from the Unstable Branch 

A boson star on the unstable branch can migrate to a new configuration of smaller 
central field density that lies on the stable branch. This behavior is analogous to the 
behavior of relativistic neutron stars • In both cases there is a large initial expansion 
of the star followed by oscillations about the final stable star configuration it will settle 
into. In the case of boson stars these oscillations damp out slowly, similar to the case of 
stars with a perfect fluid equations of state. The neutron star with perfect fluid equation 
of state settles into a configuration on the stable branch of slightly lower mass than the 
initial unstable star. The initial zero temperature star is heated due to gravitational 
binding energy being converted to internal energy by shock heating. The boson star 
on the other hand loses mass to scalar radiation and settles to a configuration of lower 
mass on the stable branch. 

The mass of a critical boson star (Model 2, Table was lowered by about 
3.5% through a radial perturbation. This large perturbation reduced the boson 
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Figure 11. The Newman-Penrose waveform ^4 is shown for a migrating critical star 
(Model 2). The initial mass of the unperturbed star was 0.633M|j/m and after a 
perturbation of the form <I>(r) — * 0.98$(r) was lowered to 0.611-Mpj/m. On top of 
this spherical perturbation a nonradial parabolic weighted perturbation proportional 
to 0.016 Y20 was applied. This was centered at about half the radius of the critical star. 
A grid size of 164 3 with resolution Aa; = Ay = Az = 0.2 was used. This waveform is 
significantly different from the nonradial quasinormal modes of a critical star. This is 
expected since the large initial perturbation moved the star significantly away from the 
original critical configuration. Subsequently it migrated to the stable branch, moving 
it further away from its initial configuration. 

field everywhere to 98% of its original value. Next, a nonradial parabolic weighted 
perturbation (see Sec. 13. lj) was superimposed on the radial one. The weighted 
perturbation was centered at a radius of R p = 3.025, which is about half the 95% 
radius of the star. This was a pure I20 perturbation with €20 = 0.016 in Eq. fTTj) . The 
data was then passed through the IVP and its evolution was tracked using the 3D code. 

The strong dynamics within the simulation causes substantial coordinate drifting. 
A dense unstable configuration with a small radius requires good resolution. This star 
migrates to a dilute stable configuration with a larger radius (about twice as large). 
The different scales place demands on both the resolution and gauge choice within the 
simulation. The gauge control was obtained by enforcing maximal slicing through a 
K-driver restoring term. When K is perturbed away from zero, this term drives it back 
exponentially in time [26J. For computational efficiency, the enforcement of maximal 
slicing is performed intermittently, with the algebraic 1+log slicing condition used on 
intermediate timesteps. While some gauge drift persists the coordinates within the 
simulation are stable. Future work will involve studying better gauge control through 
an adaptive mesh refinement technique [13*] . 

Fig. fTTI shows ^4 at a detector at r = 30. The star is dynamic and moves away from 
the initial state even in the short damping timescale of the gravitational wave signal. 
Hence, the gravitational wave does not carry the signature of the original configuration. 
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The migration process to a stable star is an extremely lengthy process, and the waveform 
shown is emitted during the very early stage of this process. Thus, the waveform does not 
show the quasinormal modes of the eventual final state stable star or the initial critical 
configuration, but rather a complex combination of the modes of a mixture of states 
that the system goes through during the migration. A good modefit of this data with 
the quasinormal modes of any single star configuration is therefore not possible. The 
lengthily and complicated precursor (both low and high frequency parts) is indicative 
of the complex dynamics of the strong perturbation and migration process. A Fourier 
analysis of this waveform in the ring down region results in a dominant frequency of 
4.3m/Mpj. This frequency is reflective of the dynamics of the star under a strong 
perturbation, which does not result in excitation of quasinormal modes. It is much 
higher than any of the quasinormal modes of a critical star (as determined by Ref. (19J). 

Since the detectors have to be far away for \l/4 to represent the gravitational 
signature of the star, a large grid had to be used to obtain the waveform. A very 
long run with such a grid would be computationally expensive. However, due to the 
strong damping of these waves, one can see the full gravitational ringing of the star 
on a short time scale (t = 60 for this simulation). On the other hand, if one wants 
to see the star continue its migration for a few radial oscillations, a longer run needs 
to be conducted. The period of one radial oscillation for this configuration is a lot 
larger than the damping time of the gravitational wave (almost four times as large). 
A two pronged approach was used. A shorter run with a large grid (164 3 of resolution 
Ax = Ay = Az = 0.2) was performed to obtain the gravitational wave. Subsequently, a 
smaller grid (96 3 of resolution Ax = Ay = Az = 0.2) was used to handle the migration. 
This was compared to a spherical perturbation run with the same grid and resolution. 

The simulation proceeds accurately for several hundred oscillations of the 
underlying scalar field, with the L2-norm of the Hamiltonian constraint on the grid 
persisting at a value of order 10~ 4 at late times. At the time at which the simulation 
is stopped, it has not reach an error prone state at which it must be terminated, but 
rather we simply exhaust our available computer time. By all available indicators, the 
simulation would continue further at this level of accuracy. The level of the errors in 
our migration study compares very favorably to other published work. 3D numerical 
simulations by Font et al. follow long-term evolutions of relativistic neutron stars 
and show errors that grow to a level of 10 -2 by the end of the simulation. The constraint 
values we report persist at least an order of magnitude below this level for hundreds of 
oscillations of the underlying scalar field, the relevant time scale for our problem. 

Fig. El shows the maximal radial metric Qi'f clS £L function of time for a migrating 
(Model 2) critical star. The initial drop in the metric indicates a migration to the 
stable branch. The star then oscillates about the new stable configuration that it will 
settle into with the radial oscillations damping out slowly. The star settles into a fairly 
constant oscillation frequency. The last four peaks are at t — 909, 1135, 1361, and 1587, 
all separated by At = 226. The metric in Fig. El also exhibits a shift upwards. We 
believe this a nonphysical effect, likely a result of gauge drift that occurs in response to 
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Figure 12. The maximum value of metric function g rr is plotted against time for 
the migration of a critical star to the stable branch. The initial mass of the star was 
0.633Mp[/m and after a perturbation of the form $(r) — » 0.98<l>(r) was lowered to 
0.611Afp[/m. The nonradial oscillations have damped out early on (approximately by 
time t = 60) and the remaining spherical oscillations proceed through seven periods 
within this simulation. The frequency is very steady by the end with the last four 
peaks at t = 909, 1135, 1361, and 1587, all separated equally by At = 226. A 96 3 grid 
of resolution Ax = Ay — Az — 0.2 was used. 

the very large initial perturbation. One can roughly estimate the final configuration to 
which the star will settle. The final star configuration of mass 0.59Mp,/m corresponds 
to a maximum metric of g rrmax = 1.18. 

Since the star expands to a configuration of larger radius the grid size had to 
accommodate this expansion. This was at the expense of the grid resolution of the 
initial denser, more compact configuration. Similar problems were seen for neutron star 
simulations. In the neutron star case a 96 3 grid was used to study the migration and the 
low grid resolution of the initial configuration caused a nonnegligible deviation of the 
average central rest mass density from the expected value [12] ■ We plan to investigate 
the use of fixed mesh refinement [43J to accommodate the multiple length scales that 
occur within such migration studies. 

In the nonradial perturbation case the high frequency oscillations, associated with 
the nonradial modes, of the peak of the metric function superimposed on the 

low frequency radial oscillation. This superposition is not clear in the figure because 
the small amplitude of nonspherical oscillation is suppressed by the larger spherical 
oscillation amplitude. The spherical oscillation amplitude is of order 0.1. This is much 
larger then the corresponding radial oscillation amplitude due to discretization errors (of 
order 0.003) from Fig. Efa) for the same timescale. While the nonspherical oscillation 
can not readily be seen, the metric exhibits the same superposition of the radial and 
nonradial frequencies discussed in Fig. |Sfa). 
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4. Conclusion 

This is the first paper to present the long time evolutions of boson stars under nonradial 
perturbations in full 3D general relativity. The existence of rapidly damping quasinormal 
modes as discussed by Ref. ^H] is confirmed. For the first time gravitational waveforms 
(both Zerilli and Newman- Penrose scalar ^4) have been fully extracted and presented 
for boson stars in 3D. We find that when a small nonradial perturbation is applied to 
a boson star, the star loses mass through gravitational radiation and has quasinormal 
modes frequencies that are much higher than its radial quasinormal mode frequency. 
The metric acquires this frequency signature. If the small nonradial perturbation in the 
scalar field is doubled the metric perturbation also doubles. For perturbations that have 
radial and nonradial components the metric exhibits both frequencies superimposed on 
one another. The gravitational wave signal is observed to damp out quickly and the 
star becomes spherical on a short timescale. 

In this paper, the gravitational waveforms have been presented for stable, 
critical and unstable boson star configurations. Our results were compared to linear 
perturbation results of Ref. ^Hj- However, the comparison could only be approximate 
because their WKB formulation leads to inaccuracy for the calculation of the lowest 
modes. Furthermore, the calculation of the quasinormal frequencies in Ref. ^5] is 
slightly dependent on the choice of the surface of the star. In principle, a boson star 
is of infinite extent with exponential damping of the scalar field within a short radius. 
This exponential damping allows one to consider a large part of the star to be virtually 
a vacuum and arbitrarily choose a surface somewhere in this region for calculational 
(perturbation theory) or numerical (finite grid) purposes. Different choices of surface 
can lead to slightly different results for the modes (changes them within a few percent). 

We have considered perturbations that are purely nonradial and are proportional to 
the spherical harmonics. The Zerilli gravitational waveforms of the i = 2, m = mode 
for small perturbations proportional to Y20 have the same frequencies as waveforms 
generated from perturbations composed of a linear combination of Y20 an d Y22 for a 
given star configuration. This shows that a given star has a specific quasinormal mode 
signature. 

It is more challenging to obtain gravitational waveforms when matter is present in 
the spacetime. Even if the system has an approximate Schwarzschild exterior region 
the dynamics of the simulation might carry matter there. The presence of matter can 
affect the Zerilli waveform which requires a region that is close to Schwarzschild. This 
effect was tested explicitly by using different types of perturbations. The first was purely 
nonradial, while the second was more general and had radial and nonradial components. 
The former did not result in scalar radiation and the Zerilli waveform was observed to 
exhibit the high frequency quasinormal mode oscillation, while the waveform in the 
latter case was noisy and seemed affected by the presence of scalar radiation. The 
Newman-Penrose waveforms behaved better in all cases. 

Nonspherical apparent horizons were observed under nonradial perturbations when 
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the horizon formed before the full emission of the gravitational waveform. Radial 
perturbations that amplify the mass of the star are seen to accelerate the collapse. 
However, the nonradial perturbations do not affect the collapse time. The geometry 
of the horizon is seen to become more spherical as the star evolves. Even for large 
nonradial perturbations the degree of asymmetry is quite small. 

Under large radial perturbations that mimic the removal of scalar field, an unstable 
branch star can migrate to the stable branch [12] . A large radial perturbation of 
an unstable star was performed for the first time in 3D and several oscillations of 
the metric were observed with the star settling to a constant oscillation frequency. 
This is a dynamic problem that involves multiple scales. Gauge control was obtained 
using a stable implementation of maximal slicing. Improvements are under study. A 
nonradial perturbation superimposed on the radial perturbation resulted in the emission 
of gravitational waves during the migration. The migration of an unstable branch star 
to the stable branch under a general perturbation (with both radial and nonradial 
components) can be significant in the formation of boson stars ^2] • 
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